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Abstract 

Excited states are extracted from lattice correlation functions using a non-uniform prior on the model pa- 
rameters. Models for both a single exponential and a sum of exponentials are considered, as well as an alternate 
model for the orthogonalization of the correlation functions. Results from an analysis of torelon and glueball 
operators indicate the Bayesian methodology compares well with the usual interpretation of effective mass tables 
produced by a variational procedure. Applications of the methodology are discussed. 
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1 Introduction 

The best means by which to extract the mass of a state 
from its lattice correlation function remains an open ques- 
tion. Common practice still relies on visual examination 
of an effective mass table for evidence of a mass plateau 
to identify the mass spectrum of a particular class of 
operator. A single correlation function will contain con- 
tributions from more than one mass eigenstate even af- 
ter application of a variational procedure, as the lightest 
state will only dominate at long time intervals. Further 
complications arise from the statistical noise of an actual 
simulation. Here, we address the question of fitting a sum 
of exponentials to the correlation functions through the 
application of Bayesian data analysis with a non-uniform 
prior on the model parameters. 

While fitting a sum of exponentials has now been 
superseded by the spectral maximum entropy method 
(MEM) [l|;[2;llli;[i], it remains in use [i;[3] when one has 
reason to believe the data may be described by a discrete 
spectrum, as to be expected for a finite lattice [1]. The 
use of least squares, or maximum likelihood analysis, to 
fit the model parameters contains an unrecognized bias 
on the magnitude of the mass. The Bayesian method- 
ology with a non-uniform prior allows one to correct for 
that bias as well as to encode additional relevant infor- 
mation. First we will discuss our choice of prior for the 
amplitude and decay constant of an exponential fitting 
function. Considering four models of exponential decay 
appropriate for lattice correlation functions, we apply the 
methodology to a collection of torelon operators, followed 
by a discussion of how one determines the most likely 
model. We next consider an alternate approach, based 
upon modelling the orthogonalization provided by the 
variational procedure, and apply that method to both 
torelon and glueball correlation functions. We close by 
summarizing our results for the spectra and discussing 
the utility of these algorithms. 

The correlation functions used herein come primar- 
ily from a 10,000 measurement run of a lattice simula- 
tion at /3 = 6 and L = 16 for SU(2) pure gauge theory 
in Z? = 2 + 1 dimensions using the Wilson action with 
spacing a. Masses are given in terms of lattice units 
throughout. Evaluation of operators occured once every 
10 compound sweeps after sufficient thermalization, up- 
dated via the Kennedy-Pendleton heat bath algorithmj9| 
augmented with a 4:1 ratio of over-relaxation sweeps [10[ 
and global gauge transformations. Further simulation 
details and particulars of the superlink method of opera- 
tor construction are found in Ref . ll| . Here we will con- 



sider operators for the torelon constructed from Polyakov 
loops and operators for the — 0+ glueball constructed 
from square boxes of diagonal superlinks. The methodol- 
ogy applies equally to general gauge groups in arbitrary 



dimension once the correlation functions are computed; 
this paper concerns itself not with the presentation of 
new results but rather with an investigation of a new 
technique for getting those results. 

2 Choice of prior 

The essential feature of Bayesian data analysis which 
takes it beyond simple least-squares fitting is the use of 
a non- uniform prior in appropriate circumstances 12|. 
Using the language of conditional probabilities [l3], we 
write "the probability of A given B under conditions J" 



prob(yl|B;/) =p{A\iB) 
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when the background information / is unchanging, and 
one states Bayes' theorem in the context of parameter 
estimation as 

(2) 
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reading "the evidence for parameters X given data D 
equals the prior for X times the likelihood of data D 
given X divided by the chance of measuring data D". 
What we call "the evidence" is often called "the pos- 
terior" , as the normalization constant affecting nei- 
ther parameter estimation nor model selection is some- 
times called "evidence"; both "prior" and "likelihood" 
have their usual meaning. The logarithm (base e) of 
Eq. ([2|) reads Le = Lp + Ll + jj^o, where the final term 
is a constant equal to — logp'^. For independent data 
D = {Dt} indexed by t with Gaussian noise a, the likeli- 
hood factors as = Y{^{2TTa1)~^/'^ eyiY>{-Rl /2), where 

Rt = [Mt{X) — Dt\/at is the normalized residual of the 
model M, so that has one term proportional to the 
measure of fit = i?f and another which is con- 
stant. With the definition of the merit function in terms 
of the model parameters, 
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(3) 



the problem is reduced(!) to one of nonlinear global opti- 
mization, with all the attendant difficulties: just because 
a solution has not been found does not mean it cannot be 
found, and just because a (local) solution is found does 
not mean it is the global one. Short of evaluating the 
merit function over the entire prior range, one must rely 
on intuition and luck to varying degrees. One's intuition, 
encoded in the form and domain of the prior functions 
, contributes to the gradient of the log evidence in 
the limit of poor data, thereby improving the chances of 
success. 



The choice of prior [ij] represents one's background 
knowledge on the likely distribution of a parameter x G 
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Figure 1: Comparison of normalized priors for various 
functions / 

[a;o,a::i] before analysis of the current set of data, and a 
uniform prior pf = = l/{xi — xq) reduces Bayes' 
theorem to a statement of proportionality between the 
evidence and the likelihood, cx p^ . A non- uniform 
prior p^ arises naturally in many contexts, often repre- 
senting a prior which is uniform over a change of variables 
a; — > F for some integrable function f{x) = dF/dx, with 
normalization pj = Ap^f{x) for Ap = J^^ f dx such 
that J^^ pj dx = 1. Besides the uniform prior, one com- 
monly encounters the Jeffreys' prior = x uniform 
over logo; and the Cauchy distribution /^^ — 1 -f a;^ 
uniform over arctanx, and we will find it useful to con- 
sider a prior we call the double Cauchy prior, /^^ = 

[1 -I- (XalxYW^ + (x/Xf,)^] for Xq < Xa < Xh < X\, that 

mirrors the form of the Cauchy prior around a central 
region. The scale parameters Xa and Xh allow one to 
introduce "soft" limits on the parameter x well within 
the "hard" limits imposed by one's evaluation range and 
make explicit the choice of units for x. These priors are 
compared in Fig. [1] 

When fitting a single exponential y[t) — 74exp(— mt) 
to data, one commonly takes the logarithm of the or- 
dinate to yield a linear model Ly{t) — La — mt with 
parameters La and m. Identifying the intercept as a lo- 
cation parameter with uniform prior indicates a Jeffreys' 
prior for the amplitude A, and the prior uniform over the 
angle tanO = m is the Cauchy distribution. Upon nor- 
malization Ly(0) = 0, only the slope m remains, whose 
best estimate Wcst from a noisy exponential with known 
decay mgig is found by minimizing the merit function 
—Le = -I- log(l + m^) -I- #. Maximum likelihood 
implies using a prior uniform on the magnitude of the 
slope m which appears in 61 as oc 1 + (tan 6*)^, clearly 
displaying a preference for a slope (or mass) of extreme 
magnitude. In Fig. [5] we compare the estimate using 
both the uniform (a) and the Cauchy prior (b) by dis- 
playing contours of the evidence, with assumed variance 




Figure 2: Contours of the evidence for the uniform 
prior (left) and the Cauchy prior (right) as the assumed 
variance decreases from at = 10~^ in (a) and (b) to 
at — 10^^ in (c) and (d). The known mass of the test 
signal is mgig ranging from 1 to 10, and the estimate 
from the evidence is most- The uniform prior does not 
resolve the mass beyond a limit imposed by the quality 
of the data, while the Cauchy prior may underestimate 
the mass for poor data 

at — 10~^, which is not the same as the nonlinear noise 
added to make the pure exponential resemble an actual 
lattice correlation function, and with at = lO"'^ in (c) 
and (d) . We see that the effect of the non- uniform prior 
VLp 7^ is to reduce the spread of the evidence be- 
yond a value of most determined by the precision of the 
data, which in practice contributes to the gradient of the 
log evidence VLe when the data has very little to say, 
VLl 0. 

3 Analysis of torelon operators 

Turning now to some real data, in Fig. [3] we display on 
a logarithmic axis the values of the normalized corre- 
lation functions Cj(t) indexed by t in units of a and 
averaged over both spatial directions of the timeslice 
for torelon operators constructed from smeared Polyakov 
loops non-contractible around the spatial lattice, noting 
that they do not tend to zero at the largest time sepa- 
ration. Our count of operators equals 6 indexed by j, 
representing 3 step sizes and 2 smearing levels. We de- 
note by C'^ the self-correlations of the original basis of 
operators and by the correlation functions produced 
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Figure 3: Positive values of the normalized self- 
correlation functions for torelon operators displayed 
against time t and index j, where the symbol type in- 
dicates the timeslice separation. The (t) in (c) and 
(d) have been orthogonalized by a variational procedure, 
while Cj^{t) in (a) and (b) are the correlations of the 
original basis. The ground state is observed from the 
regularity in spacing for j = 1 in (d) 



by a variational procedure to enhance their orthogonal- 
ity 15; 3 121- Specifically, we take the modal matrix E 
of real eigenvectors found by diagonalizing the product 
of the first two cross-correlation matrices Xt = Xij (t) = 
{0*{t + T)Oj{T))T for zero-momentum timeslice opera- 
tors 0{t) such that 



X(, ^XiE 



ED 



(4) 



where D is a diagonal matrix arranged by increasing 
eigenvalue. Symmetrizing the cross-correlations on i,j 
is equivalent to averaging them over both temporal ori- 
entations ±t, so that Xt -l-X^ = Xf -|-Xi_t, and we will 
work with the symmetrized data with temporal range 
t £ [Q,L/2]. The orthogonal basis OE is formed from 
the original operators O = Otj ~ Oj{t) so that Vf — 
E-^XjE, and the variational self-correlators are extracted 
Cj{t) = 5,,V,,{t). ^ 

Common practice is to form the effective mass ta- 
ble from the correlation functions, defined by amj{t) = 
\og[Cj{t — 1)/Cj{t)] for t > 1, with errors given by jack- 
knife analysis. In Table [T] we display the effective mass 
table for both the original and variational torelon op- 
erators, with negative and imaginary values zeroed out 
and neglecting the error analysis. Clear evidence for a 



Table 1: Effective mass table for original am^ and vari- 
ational amj torelon operators 



t 


1 


2 


3 


4 


5 


6 


7 


8 




1.34 


1.10 


1.06 


0.98 


0.85 


1.40 


0.00 


0.09 




1.28 


1.09 


1.06 


1.01 


0.84 


1.58 


0.00 


0.00 
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1.18 


1.06 


1.07 


1.05 


0.92 


2.23 


0.00 


0.00 


a 


1.16 


1.05 


1.07 


1.05 


0.94 


2.41 


0.00 


0.00 




1.18 


1.06 


1.07 


1.04 


0.92 


2.24 


0.00 


0.00 




1.07 


1.01 


1.05 


1.03 


1.12 


0.00 


0.00 


0.00 




2.53 


1.87 


1.49 


0.00 


0.56 


0.31 


0.00 


0.00 




3.52 


4.43 


0.00 


2.11 


0.00 


0.60 


0.00 


0.00 




4.26 


3.16 


0.67 


0.00 


0.00 


0.00 


1.34 


0.00 




4.63 


3.14 


0.83 


0.00 


0.00 


0.00 


2.36 


0.00 




5.19 


2.41 


0.00 


2.66 


0.00 


0.00 


0.00 


0.03 



mass plateau at the ground state identified by the vari- 
ational procedure is seen for the correlation functions 
of the original basis, but the remainder of the data is 
hard to interpret — how should one identify and extract 
the mass of the excited states? In other words, how is 
the information obtained for the remaining states (which 
here appear to be multi-torelon excitations) by the vari- 
ational procedure? To address these questions, we con- 
sider fitting a sum of exponentials with free parameters 
for the amplitudes and decay constants, a notoriously 
hard problem [3| whose difficulties we hope to mitigate 
through the use of non-uniform priors. 

3.1 Models for exponential decay 

The first of the models we consider, conveniently indexed 
by their number of parameters, is given by a single ex- 
ponential with a free decay constant. 



Ml : C{t) = exp(-TOii) , 



(5) 



which is driven primarily by the first non-constrained 
value C(l). As an alternative, we consider a model which 
utilizes a constant to represent the statistical noise of the 
simulation, 



M2 : C{t) = Ai exp(-TOit) + (1 - ^i) 



(6) 



chosing a double Cauchy prior for mi in both models 
with hard limits of 0.1 and 6 and soft limits of 1 and 4 
in lattice units of mass. (With such a tight range, the 
shape of the double Cauchy prior approaches that of an 
offset Gaussian on logarithmic axes.) For the amplitude 
Ai, the Jeffreys' prior is chosen over range [0.8, 1.1]. De- 
tailed observation of the correlation functions indicates 
that fluctuations have become dominant several sites be- 
fore the midpoint of the lattice is reached, so we restrict 
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Table 2: Single exponential model parameters for varia- 
tional basis torelon operators (t) 



Ml 


Ah 


ami 






am\ 




1.071(4) 


-11.2 


0.997(2) 


1.077(6) 


-11.1 


2.534(22) 


-9.9 


0.997(1) 


2.572(26) 


-10.7 


3.524(60) 


-14.3 


1.000(1) 


3.514(67) 


-9.9 


4.257(125) 


-10.7 


1.000(1) 


4.236(137) 


-9.8 


4.625(180) 


-7.6 


0.999(1) 


4.684(212) 


-10.3 


5.172(310) 


-8.3 


1.001(1) 


5.058(312) 


-11.0 



the fitting window to a range i^t G [0,5] and neglect to 
replace the exponential with a hyperbolic cosine repre- 
senting correlations the long way around a finite lattice, 
justified when mL 1. (The variance of the correla- 
tors at < = is calculated before the normalization is 
applied.) In Tabled] we give the results for these two sin- 
gle exponential models applied to the variational basis 
correlations Cj{t) with the standard error indicated in 
parentheses and Vio^Z^b = log^o IVL^;]. 

To accommodate the covariance in t of the functions 
Cj{t), one generalizes the measure of fit by 



{M,-C,fW,iM,-C,), 



(7) 



where Mj and Cj are column vectors indexed by t and 
W = is the inverse of the variance matrix evalu- 
ated [8| from the iV„ configurations indexed by n, 

s ^ = (7v„ - i)-i((cr - {cix){c: - {c:x))^ . 

Note that it is the matrix inverse jl2| which appears in 
Eqn. ([T]) , a distinction lost in the more common notation 
of the general measure of fit as a sum over indices. 

Considering now a sum of two exponentials, the third 
model is just the simple sum without a noise floor, 

M3 : C(i) = ^iexp(-mit) + (l-Ai)exp(-m2i) , (9) 

and the fourth has two free amplitudes with a noise floor, 

C{t) — Ai exp(— TOit)-l-A2 exp 



m2t)+{l-Ai-A2) 
(10) 

We maintain the double Cauchy prior on the masses and 
the Jeffreys' prior on the amplitudes, this time with range 
[0.1, 1]. These models are applied to the original basis of 
correlation functions C^{t). Their solution parameters 
are displayed in Table ^ We see that the original self- 
correlators do contain evidence of at least one excitation; 
however, the resolution is poor in terms of the standard 
error. When compared graphically in Fig. |31 where for 
the sum of exponentials models the secondary states are 




Figure 4: Comparison of the spectra of the torelon oper- 
ators indexed by state j for the various models 



marked by open circles while the dominant states are 
filled circles, we find a consistent estimate for the ground 
state, whereas the estimate for the excited states depends 
upon the use of the noise floor constant when evaluated 
from Cf{t). 

3.2 Model selection 

How does one compare the quality of fit between the 
models? The Bayesian formalism addresses model se- 
lection by considering the ratio of the evidence for each 
model with no prior preference 
to the likelihood ratio 



p , thus reducing 



Pd 
Pd 



Pa 
Pb 



(11) 



whose factor for each model is the "integrated probability 
bump" over the model parameters X, 



pZ 



D.X, ^ 

Pm 



PmP°x.aM 



(12) 



(normalized by the chance of D which factors out of the 
ratio). Our choice of nomenclature [c/. Eqns. ^ and 
(fTT|) ] identifies Eqn. (fT2|) as "the likelihood of D given 
M" (not "the chance of D given M" ) which marginalizes 
into a product of the prior for X and the likelihood for 
X. In other words, the evidence for the model is the (un- 
normalized) integral of the evidence for its parameters. 
Introducing a model subscript M to Eqn. one sees 
that is the normalizing constant such that / 
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Table 3: Double exponential model parameters for original basis torelon operators Cf {t) 



Ms 


Mi 


Ai 








Ai 


A2 








0.654(137) 


1.020(81) 


2.580(769) 


-11.5 


0.751(87) 


0.248(88) 


1.089(63) 


3.487(1670) 


-10.5 


0.732(128) 


1.033(70) 


2.655(993) 


-10.9 


0.804(65) 


0.195(66) 


1.082(47) 


3.598(1782) 


-9.9 


0.614(162) 


0.993(95) 


2.316(600) 


-11.0 


0.754(104) 


0.244(105) 


1.087(71) 


3.257(1579) 


-10.9 


0.851(61) 


1.042(34) 


3.242(1653) 


-11.7 


0.869(37) 


0.131(38) 


1.054(30) 


3.846(1918) 


-10.1 


0.868(55) 


1.039(31) 


3.268(1714) 


-11.5 


0.883(34) 


0.117(35) 


1.047(28) 


3.840(1929) 


-10.8 


0.855(63) 


1.041(35) 


3.182(1644) 


-12.0 


0.874(37) 


0.125(38) 


1.053(29) 


3.824(1919) 


-10.5 



Table 4: Negative log evidence of the models for torelon 
operators 



model 




equals unity. An unfortunate confusion of nomenclature 
arises because p^.j appears both in the position of chance 
in Eqn. ([2]) and in the position of likelihood in Eqn. ([TT|) . 

Under the quadratic approximation, generally accept- 
able when the evidence is not severely truncated by the 
prior range, one can evaluate the integral analytically to 
write the negative logarithm of the likelihood as 



L 



1 



k 



log/-i- Vlog A-i^W)^ 



fc 



(13) 

for X indexed by k and {hk] the eigenvalues of the in- 
verse variance Jlfc ^fe ~ detS^^, where the first two 
terms are the value of the merit function evaluated at its 
minimum and the remainder comprise the Occam factor 
accounting for the ratio of the width of the evidence S 
to the prior range {A^}. An additional parameter must 
provide not just a better fit but a significantly better fit 
in order for its plausibility to increase. With several mod- 
els to choose from, the one with the lowest value of —L^.[ 
is deemed the most plausible, with the relative proba- 
bility given by the exponential of the difference between 
the (negative) log evidence for each. (Some investigators 
take the further step of normalizing X)mPm ~ ^ which 
we neglect.) In Table S] we display for our four 

models, noting that, except for one case, the models Mi 
and M3 without the noise constant are preferred over 
those with its inclusion. 




Figure 5: Effective masses from singular value decom- 
position for torelon (a) and 0"*" glueball operators (b) 
indexed by time-slice separation t, where only positive 
eigenvalues are marked. The symbol type indicates the 
ordering of the states j , and the estimate for the ground 
state o is observed from the lowest horizontal values 

4 A better approach 

The methodology of the previous section, while produc- 
ing a credible analysis of the torelon data, did not fare 
well when faced with glueball correlation functions. Con- 
sequently, we have investigated an alternate approach 
which takes into account the orthogonality of the mass 
eigenstates. The focus of our investigation is a compari- 
son of the two means by which one may estimate param- 
eters from a set of data, either by inverting the data or 
by inference from the data [HI; [3] • There are a variety 
of ways to invert a set of correlation functions to produce 
effective masses but let us consider one based upon 
the variational method mentioned previously. Rather 
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than using the eigenvectors of Eqn. (jlj, one can apply 
the diagonahzation in succession to the original basis of 
cross-correlators for t > 0, 



Xq ^XjEt — EtDt 



(14) 



and then define effective masses in terms of the eigenval- 
ues amj^{t) = — t^^ log Djj{t). One finds that the same 
eigenvalues are produced by a two step procedure where 
one first orthogonalizes the basis according to the zero 
separation correlators, 



Yt = (AoDo '^YXtAoDo 



(15) 



where Xg = AoDqAJ^ so that Yq = I, and then in- 
verts the data Yt = EjDfE^ to get the effective masses 
amj (t) = -~t~^ log Djj{t). As Yt is symmetric, one may 
as well use singular value decomposition Yt = ASB-^, 
where the columns of A and B may differ only by a sign 
for negative eigenvalues Djj — —Sjj, to keep the values 
ordered by magnitude for display. We see in Fig. [5] how 
the mass eigenstates peel away at small t from the band 
of values we interpret as simulation noise. The essen- 
tial difference between Eqns. and (|15l) is that the 
orthogonalization of the original basis Xt is done asym- 
metrically in ((T4)) . coming in from the left in the guise 
of Xq ^, whereas in (|15p all Yt are symmetric. It is from 
data Yt that we shall infer values for the mass and am- 
plitude parameters. 



4.1 Modelling the variational procedure 

To model the orthogonalization provided by the varia- 
tional procedure, we recognize that the step Xt — Yf 
produces an orthogonal basis of operators which are not 
necessarily mass eigenstates. The estimates for the mass 
and amplitude parameters should have no temporal de- 
pendence, and so we write the model as 



Yt = ADtA^ 



(16) 



where A-^A = I and Dt is the decay matrix Djj{t) — 
exp{—tmj). In principle, the model can handle extract- 
ing a subset of eigenstates with the inclusion of some 
identity matrices. 



Yt=A(Dt-Ife)A^ + I<t 



(17) 



for k < d and A of d rows and k columns, but we did not 
have much success with its convergence and will focus 
on the case k = d for a subset of the entire data, Yij{t) 
for i,j < d, ordered by decreasing value of their t = 1 
self-correlator. 

While there are many options for representing the 
orthogonal amplitude matrix A, we choose to use the 
minimal number of degrees of freedom provided by the 



composite parametrization presented by Spengler et 
al. Summarized for the real case here, one starts with a 
vector O of d{d — l)/2 angles with domain 8; G [— tt, tt) 
mapped into = 6„m using / = n — d + m{2d — m — l)/2 
for 1 < m < n < d. Using bra-ket matrix notation, the 
generators are defined by Amn = ~i\m)(n\ + i\n){m\. 
The amplitude matrix is then given by the ordered prod- 
uct 

d-i r d 

A ]J Yl exp(ie„i„A™„) , (18) 

m— 1 [_n—m-\-l 

and a uniform prior is selected for the angles Qmn- In 
order to keep the mass estimates distinct, we parametrize 
their energy gaps above the vacuum Eq = such that 
nij = J2i<j fo'' j > 0. In terms of the energy gaps, the 
decay matrix reads Dt = '3xp(— t _Bi)|j)(j|. The 
double Cauchy prior is chosen for the gap parameters, 
the first with hard limits of 0.1 and 4 and soft limits of 
0.5 and 1 (1.5 for the glueball operators), and the rest 
with limit pairs of [0.01, 2] and [0.1, 1]. 

4.2 Fitting the data 

This parametrization of the orthogonalized cross-correlators 
is unique only up to some leftover real phases, Yt = 
ADtA^ = BDtB^ where B = AI for I = J2k Mk){k\- 
Our optimizer did not like the restriction of the param- 
eter domain to [0,7r), so we address the redundancy of 
the model by including a multi-modal factor of 2^^ on the 
integrated probability bump. The measure of fit is now 
a matrix quantity. 



which we reduce to a scalar using the 1-norm 



llx: 



d d 

EE 

1 = 1 J = l 



(19) 



(20) 



This formula is identical to that for the variance of a sum 
of quantities with covariance [l^l which we will need to 
give error bars to our estimates for the eigenstate masses 

The best fitting parameters for our 10,000 measure- 
ment (10k) data are displayed in Table [5] (with the stan- 
dard error below the parameter value) for both torelon 
and 0''" glueball operators. We have restricted these Yf 
to d = 4 to extract the lightest states, and the entire 
range of t is included in the fit. We also show the pa- 
rameters for an independent run of 1,000 measurements 
(Ik) using a similar but updated code. (The Ik run av- 
eraged correlators over the D — 3 independent temporal 
orientation thus is more like 2,500 measurements of the 
10k run which did not.) With only one model per set 
of data, there is no opportunity for model selection, but 
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Table 5: Best fit parameters for torelon and 0+ glueball operators with the standard error below the parameter 
value 



Nn 


o 


El 


E2 


Es 


E4 


61 


62 


63 


64 


65 


66 


VlO-^B 


-Le 


O 


tor 
a 


1.083 
0.004 


1.440 
0.018 


1.659 
0.094 


0.509 
0.176 
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Table 6: Mass and amplitude parameters for torelon and 0+ glueball operators 










10k 








Ik 








3 


1 


2 


3 


4 


1 


2 


3 


4 






1.083 


2.523 


4.182 


4.691 


1.141 


2.779 


3.500 


4.471 




a 


0.004 


0.018 


0.096 


0.165 


0.008 


0.042 


0.095 


0.234 


5-1 




0.915 


-0.400 - 


-0.046 


0.018 


0.891 


-0.448 


-0.038 


0.059 




A 


0.331 


0.679 


0.645 


-0.113 


0.342 


0.568 


0.203 


-0.721 






-0.147 


-0.415 


0.368 


-0.819 


-0.277 


-0.658 


0.474 


-0.515 






-0.176 


-0.455 


0.668 


0.562 


0.112 


0.210 


0.856 


0.459 






1.271 


1.749 


2.347 


2.815 


1.278 


1.841 


2.440 


2.843 




a 


0.005 


0.007 


0.015 


0.027 


0.013 


0.022 


0.040 


0.066 




0.832 


0.512 - 


-0.208 


0.054 


0.825 


0.478 


0.299 


-0.042 




A 


-0.494 


0.520 - 


-0.607 


0.342 


-0.522 


0.445 


0.699 


-0.202 






0.248 


-0.622 - 


-0.372 


0.643 


0.209 


-0.671 


0.417 


-0.576 






0.052 


-0.285 - 


-0.671 


-0.683 


0.062 


-0.350 


0.498 


0.791 



8 




0.5 \ 
0.1 I r 



O o 



-0.1 
0.02 



-0.02 
0.01 



-0.01 



2 4 6 
t 




0.01 



-0.01 
0.01 



-0.01 



• O Q O 6 -0.1 

(0) 1 



(d) 



(h) 



(I) 



■ O O 1 > 



0.5 





2 4 6 
t 



(P) 





2 4 6 

t 



Figure 7: Comparison of Yt = Yij{t) for the model (solid) and data (circled) for the 0+ glueball operators using 
10k data 
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Table 7: Best fit parameters using Ik data for 0+ glueball 
operators including the vacuum contribution with the 
standard error befow the parameter value 
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1 



0.000 
0.000 



1.271 
0.008 



0.582 
0.023 



0.593 
0.045 



0.426 
0.075 



1 



-0.036 
0.000 



0.029 
0.000 



-0.022 
0.000 



-0.015 
0.000 



0.631 
0.015 



6 



10 



a 



-0.293 -0.147 0.855 0.416 -0.572 
0.013 0.010 0.037 0.028 0.087 



we include the value of the negative log evidence in the 
table for completeness. 

The parameters of our model must be manipulated 
to produce estimates for the quantities in which we are 
interested, namely the masses and amplitudes of the en- 
ergy eigenstates. As the mass nij is the sum of the en- 
ergy gaps Ei<j, its variance is given by the 1-norm of 
their covariance, (t,^„ = As there must be a high 

degree of correlation in A to enforce its orthogonality, 
and as the formula for — ^ A is rather complicated, we 
have not computed standard errors for the entries of A. 
In Table [H] we give the mass and amplitude parameters 
for our 10k and Ik runs. These amplitudes, one recalls, 
are the overlaps with the orthogonalized Yt data which 
themselves are a linear combination of the original basis 
Xt . We note the appearance of two negative real phases 
in Abox between the 10k and Ik data for the states la- 
belled 3 and 4 which explains why the final value of Obox 
differs in Tabled 

Finally, there is the graphical inspection of the quality 
of fit. In Figs. [5] and [7] we compare the of the model 
and the data for the torelon and 0+ glueball operators. 
For either operator the model and data are in visibly 
close agreement, except for Y^^^{t) which appears to be 
mostly noise. The significant variation in the model lies 
primarily within the first quarter lattice extent, which 
here is spanned by only 5 data values — we simply are 
not working at sufficient temporal resolution to measure 
the cross-correlations where they vary most. We also see 
that (for niL ^1) the relevant part of the data for fitting 
the model parameters lies well away from the midpoint 
of the lattice thus should not be much affected by the 
difference between exp and cosh. As a rule of thumb, if 
one's results depend upon the use of cosh, then one is 
not working with a large enough lattice. 



Table 8: Mass and amplitude parameters using Ik data 
for 0^ glueball operators including the vacuum contribu- 
tion 
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0.000 
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a 


0.000 


0.008 


0.022 


0.041 
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0.999 
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-0.005 


-0.001 




0.036 


0.763 


0.562 


0.313 


-0.054 


A 


-0.029 


-0.559 


0.332 


0.735 


-0.190 




0.022 


0.286 


-0.644 


0.349 


-0.618 




0.015 


0.147 


-0.400 


0.489 


0.761 



4.3 Fitting the vacuum 

Our final application of the composite parametrization 
of the orthogonalized cross-correlators addresses whether 
one needs to work with vacuum subtracted operators 
for the O"*" glueball or not. We note that the need 
not be normalized, as Eqn. (jl5p produces a normalized 
Yt by construction. Neither must one take its vacuum 
subtracted value, as the vacuum is simply the ground 
state of operators with trivial quantum numbers. To fit 
the vacuum contribution to Yt without vacuum subtrac- 
tion, we simply prepend the first gap parameter with 
one whose prior is uniform over ±0.1 and treat the first 
nonzero energy state as an excitation. With this analysis 
we get for the Ik run of data the best fitting parame- 
ters shown in Table [71 The vacuum energy evaluated to 
Eq = 0.32(35) X 10"^, and this fit's quality is given by 
Vio^B = -7.4 and -Le = 209. From Table [S] we see 
that the mass estimates for states above the vacuum are 
consistent with those in Table[6]using vacuum subtracted 
correlators. 



5 Summary and conclusions 

We summarize our results for the mass eigenstates of the 
torelon and 0+ glueball on a 16^ lattice at /3 — 6 us- 
ing the composite parametrization of the orthogonalized 
cross-correlators in Fig. |S1 While there is a little bit of 
shifting of the estimates for the excited torelon states 
between the 10k and Ik data, the estimates for the 0+ 
glueball are consistent for the long and short measure- 
ment runs. The expectation of a discrete spectrum is 
observed for these lightest few energy eigenstates. While 
there is an improvement in the error bar width for the 
10k data, the estimate from the Ik data is already fairly 
precise. (One of course should not confuse precision with 
accuracy.) By obviating the need to normalize and vac- 
uum subtract one's timeslice correlators, there should be 
a modest decrease in one's execution time (not investi- 
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Figure 8: Mass spectra for torelon (a)-(b) and 0+ glueball 
operators (c)-(d) using the composite parametrization 
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